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ABSTRACT 


A  mixed  system  of  parabolic  and  elliptic  partial  differential  equations 


is  used  to  describe  the  carrier  transport  and  potential  distribution  in  semi- 
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conductor  devices  such  as  MOSFET's,  diodes  -evt-rcrs  A  singular  perturbation 

* 


analysis  of  the  corresponding  initial  boundary  value  problem  is  carried  out. 


Asymptotic  expansions  of  the  solution  in  powers  of  the  minimal  Debeye  length 


are  given.  Based  on  this  analysis  a  finite  difference  method  for  the 


numerical  solution  of  these  problems  is  developed.  Here  problems  arise  due  to 


different  time  scales  which  are  intrinsically  present  in  the  analytical 


problem.  These  different  time  scales  do  not  occur  in  the  physical  solutions 


because  of  special  (equilibrium-)  initial  conditions.  Nevertheless  they  cause 


severs  stability  problems  for  finite  difference  methods.  An  unconditionally 


stable  scheme  is  developed  which  minimizes  computational  effort.  Numerical 


experiments  on  a  test  problem  in  one  space  dimension  are  presented. 
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SIGNIFICANCE  AND  EXPLANATION 


In  this  paper  we  are  concerned  with  a  one  dimensional  transient  model  for 
a  simple  p  -  n  junction  (i.e.  a  diode).  The  model  consists  of  an  initial 
boundary  value  problem  for  a  singularly  perturbed  nonlinear  system  of  elliptic 
and  parabolic  P.D.E.'s  in  one  space  and  one  time  dimension.  One  of  the  main 
purposes  of  this  model  is  to  give  information  about  the  "rise-time"  of  the 
p  -  n  junction  (the  time  required  to  reach  a  steady  state).  We  carry  out  our 
asymptotic  analysis  of  the  problem  for  large  doping  concentrations.  (This  is 
the  case  for  modern  devices.)  We  derive  asymptotic  expansions  of  the  solution 
in  powers  of  a  singular  perturbation  parameter  which  is  proportional  to  the 
minimal  Debeye-length  of  the  device.  It  turns  out  that  the  solutions  vary  on 
a  time  scale  which  is  proportional  to  the  average  lifetime  of  electrons  and 
holes.  Therefore  the  risetimes  will  be  of  the  same  order  of  magnitude  as  this 
lifetime.  In  addition  to  this  timescale  there  is  a  second  "first"  timescale 
Intrinsically  present  in  the  problem.  This  first  timescale  -  although  never 
observed  physically  because  of  special  (equilibrium-)  initial  conditions  - 
causes  severe  stability  problems  in  the  numerical  solution  of  the  problem. 
Based  on  our  asymptotic  analysis  we  derive  an  unconditionally  stable  finite 
difference  method  for  the  solution  of  the  p  -  n  junction  equations. 


NUMERICAL  METHODS  FOR 
TRANSIENT  SEMICONDUCTOR  DEVICE  MODELLING 


Christian  Ringhofer 

1*  Introduction 

In  this  paper  we  present  an  asymptotic  analysis  of  a  singularly  perturbed 
mixed  system  of  parabolic  and  elliptic  differential  equations  modelling  a 
p  -  n  junction.  The  physical  situation  is  as  follows:  A  semiconductor  (for 
instance  silicone)  is  doped  with  acceptor  atoms  in  the  left  hand  side  and 
donor  atoms  in  the  right  hand  side  and  a  bias  U  *  UA  -  Uc  is  applied  to  the 
contacts. 


cathode  contact 


Uc  t  Applied 

Cathode  potential 


•3 

The  device  is  assumed  to  have  a  characteristic  length  2JL  (~  0.5*10  cm)  and 
the  junction  is  at  x  *  0,  The  device  is  forward  biased  for  U  >  0  and 
reverse  biased  for  U  <  0.  The  physics  of  a  p  -  n  junction  is  explained  in 
Sze  [1969],  Ashcroft  et  al  [1976]  and  R.  A.  Smith  [1978].  The  equations 
describing  the  electrostatic  potential  and  the  carrier  and  current  densities 
were  first  given  by  Van  Roosebroeck  in  [1950].  In  the  case  of  one  space 


dimension  they 

consist  of 

(1.1) 

♦xx  -  f  (n  -  P  -  C) 

(Poisson's  equation) 

(1.2) 

Jn  "  *(Dnnx  -  wnn*x) 

electron  current  relation 

(1.3) 

Jp  "  -*(Vx  +  V**) 

hole  current  relation 

Sponsored  by  the  United  States  Army  under  Contract  No.  DAAG29-80-C-0041 


(1.4) 


continuity  equations 


(1.4)  n.  -  ± 

t  q 

<?n  ♦  K  continuity  equations 

X 

d.5)  Pt  -  - 

—  J  +  R  for  electrons  and  holes. 

q  px 

The  system  (1.1)  to  (1.5)  is  subject  to  the  boundary  conditions 

(1.6) 

-  UT  tn(p,.t#t|)  *  0A(t) 

(1.7) 

♦<*•*>  -  0T  »n(nU‘t))  ♦  0c(t) 

(1.8) 

n(±t,t)  p(±l,t)  -  n* 

(1.9) 

n(±Jt,t)  -  p(±t,t)  -  C(±t)  -  0 

and  tha  initial  conditions 

(1.10) 

♦(x,0)  ■  tj(x) 

(1.11) 

n(x,0)  *  nj(x) 

(1.12) 

p(x#0)  «  Pj(x) 

The  dependent  variables  (with  units)  are 

Table  1 

♦ 

electric  potential  (V) 

*x 

electric  field  (Vcsf1) 

n 

electron  density  (cm  ) 

P 

hole  density  (cm“^) 

Jn 

electron  current  density  (Acm-^ 

JP 

hole  current  density  (Acm  '  * 

The  parameters  q#  t,  un»  Pp*  Dn#  Dp#  Up  in  (1.1)  -  (1.5)  have  the 
following  Meaning  and  approxinate  values  at  T»  300K  (room  tenqperature ) 


Table  2 


q 

elementary  charge 

10'19As 

e 

permittivity  constant 

10"12As/Vcm: 

wn 

electron  mobility 

103cm2/(Va) 

MP 

hole  mobility 

103cm2/(Vs) 

ni 

intrinsic  density 

°n  D 

1010c»-3 

uT 

»  —  "  — ^  thermal  voltage 
u  u 

0.025V 

C  is  a  given  function  of  x  and  models  the  doping  profile,  i.e. 
(1.13)  C(x)  -  N*(x)  "  Vx)  (cm"3) 


where  N*  and  N~  are  the  densities  of  electrically  active  donor  and 
acceptor  atoms  in  the  doped  semiconductor.  C(x)  is  negative  for  x  <  0  (on 
the  p-side )  and  positive  for  x  >  0  (on  the  n-side ) .  We  further  assume 
that  C  is  an  odd  function  i.e. 

(1.14)  C(x)  *  -(C-x)  holds. 

(Although  this  does  not  appear  to  be  any  special  case  from  the  physical  point 
of  view  it  will  simplify  the  analysis  considerably.)  We  assume  that  C(x) 
has  a  jump-discontinuity  at  x  *  9*  R  in  (1.4),  (1.S)  denotes  the 
generation-recombination  term  and  describes  the  rate  at  which  electron  hole 
carrier  pairs  are  generated  or  recombine  (vanish).  For  our  purpose  it 
suffices  to  consider  the  Shock ley-Read-Ha 11  recombination  term 


np  -  n. 


(1.15) 


FSM)(n,p> 


t  (n+n, ) 
P  i 


Tn(P+ni> 


,  ~3  “1% 
(cm  s  ) 


T  and  t  denote  the  average  lifetimes  of  electron  and  holes e  For 
n  p 


T  “  10~6(s).  Other  mote  realistic  recombination 
rates  are  given  in  Langer  et  al  [1981]  and  Schfltz  et  al  [1981]#  In  general  it 


simplicity  we  set  x^  *  x^ 
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can  be  said  that  the  choice  of  the  recombination  rate  is  not  very  important 

for  the  behavior  of  solutions  of  the  time  dependent  problem  (whereas  it  is  in 

the  steady  state  case  -  see  Markowich  and  Ringhofer  [1983],  Markowich  et  al 

[1983a] ,  Mock  [1983]  and  Schtttz  et  al  [1982]). 

A  more  restrictive  assumption  is  =  const.  In  general  yn  and 

U  can  be  modelled  as  functions  of  the  electric  field  \b  .  This  might 

P  x 

influence  the  behavior  of  the  solutions  considerably  (see  c.f.  Mock  [1983]). 

For  the  rest  of  this  paper  we  assume  that  the  initial  conditions  (1.10)  - 
(1.12)  are  cooq^atible  with  the  boundary  conditions  and  the  differential 
equation.  So  we  assume  that 


(1.16) 


(1.17) 


(1.18) 


(1.19) 


(1.20) 


holds. 


♦  =  -^  ( n  —  p  —  C) 

wi  e  I  pi 

XX 

n^ti)  •  Pjdi)  *  n^ 
n jtti)  -  px(±*)  -  C(±i)  *  0 

*  V°> 

V‘>  -  °T  +  “o'01 


After  an  appropriate  scaling  we  carry  out  an  asymptotic  analysis  of  the 
problem  (1.1)  -  (1.12)  where  the  (dimensionless)  quantity 


(1.21) 


EO„  1/2 

T  x 


X  -( - * - 5) 

q  max | C ( x ) | t 
xet-*,t] 


acts  as  a  perturbation  parameter. 

We  investigate  the  behavior  of  the  solutions  as  X  ♦  0  which  corresponds 
to  very  high  doping  (which  is  the  case  for  modern  devices).  It  turns  out  that 
there  are  two  different  timescales  intrinsically  present  in  the  equations 


'-MMSK' 


/  / 


(1.1)  -  (1.5).  However  if  we  choose  the  initial  conditions  ,  nr  and  pj 
in  (1.10)  -  (1.12)  to  be  the  solutions  of  the  steady  state  problem 

(1.22)  *xx  "  f  (n  '  p  “  C) 

d.23)  Jn  -  q(Dnnjt  -  unn*x) 

(1.24)  Jp  -  -q(Dppx  +  ypp*x) 

(1.25)  J  -  R  ,  J  -  -R 

n  P 

X  X 

the  solutions  will  vary  only  on  the  "slow9*  time  scale  (which  is  of  the  same 
order  of  magnitude  as  T  the  average  lifetime  of  electrons  and  holes). 
Although  the  other  (Nfastn)  time  scale  is  not  present  in  this  case,  it  causes 
severe  stability  problems  for  numerical  methods  to  solve  (1.1)  -  (1.12).  We 
derive  a  finite  difference  method  for  problem  (1.1)  -  (1.12)  which  is 
unconditionally  stable  and  of  second  order  in  time.  At  each  time  step  the 
solution  of  two  coupled  nonlinear  differential  equations  is  required.  The 
generalization  of  this  method  to  nonsymmetric  p  -  n  junctions  (i.e.  C(x) 
not  odd)  and  to  two  or  three  space  dimensions  is  straightforward. 

This  paper  is  organized  as  follows.  In  Chapter  2  we  scale  problem  (1.1) 

-  (1.12)  and  reformulate  it  as  a  singular  perturbation  problem.  In  Chapter  3 
we  review  the  asymptotic  analysis  for  the  steady  state  problem  (1.22)  -  (1.25) 
given  in  Markowich  and  Ringhofer  [1982]  and  Markowich  et  al  [1982].  In 
Chapter  4  we  derive  the  asymptotic  expansion  of  the  solution  of  (1.1)  - 
(1.12).  In  Chapter  5  we  construct  the  finite  difference  method.  In  Chapter  6 
we  present  some  numerical  results  on  a  test  problem. 


2.  Reformulation  as  a  singular  perturbation  problem 

In  this  chapter  we  scale  the  system  (1.1)  -  (1.5)  appropriately  and 
transform  it  into  a  singular  perturbation  problem.  Essentially  we  use  the 
scaling  given  in  Markowich  and  Ringhofer  [1982]  and  Markowich  et  al  [1983a]. 

He  scale  the  independent  variables  x  and  t  by  the  characteristic  length  of 
the  device  and  by  the  average  lifetime  of  electrons  and  holes. 

«.’>  XB-f  ,  ts-f  . 

For  the  dependent  variables  n#  p,  Jn  and  Jp  we  employ  the  following 
scaling: 

(2.2)  ♦  -  ,  n  -  2  ,  p_  =  £ 


•  ".-5  •  *.-£ 

T  C  C 


(2.3) 


(2.4) 


J  «  — —  J  ,  J  — —  J 

n  —  _  n  p  —  _  p 

8  CqD  CqD  * 

n  p 


C  !■  max  |C(x)| 


In  this  scaling  equations  (1.1)  -  (1.5)  assume  the  form 


(2.5) 


2  i 

X  — r  ♦  *  n  -  p  -  C 
*  2  T  s  *s  s 
3x 

8 


(2.6) 

9  n  3  , 

8t.  .  ,2  JT  n. 

(2.7) 

p  - - E  j 

at  PS  .2  ST  p 

8  l  S  *S 

(2.8) 

J  *  t— -  n  —  n  i  * 1  (i 

n  3x  8  s  5x  Ti 

ss  s 

,  a  a 

(2.9) 

p  8x  Ps  Ps  3x 

*s  s  s 

eu.  1/2 

(2.10) 

»  -  •  e.. 

£  q  C 

(2.11) 

%  "  (nSP8  ■  °2>(n8  +  P. 

Ths  boundary  conditions  (1.6)  -  (1.9)  transforms  to 


♦_,<t>  -  J-  V*»  *  (n[c.-M  wc<-,,^  ^  , 


♦1,t,-l-Vt,,„(cl,>  -*») 


a  :*  ni/C 


and  the  initial  conditions 

(2.24)  +(x,0)  -  ♦*<*),  n(x,0)  »  n^x),  p(x,0)  *  p*(x) 

(2.19)  -  (2.24)  represents  a  singularly  perturbed  system  consisting  of  an 
elliptic  equation  (2.19)  and  two  parabolic  conservation  laws  (2.20)  - 
(2.21).  The  compatibility  conditions  (1.16)  -  (1.20)  now  read 

(2.25)  ^(-l)  -  ♦_■,«>),  ♦I(1)  -  f1(0) 

(2.26)  n1  -  j  (C  +  /C2  +  4o2),  p1  =  j  ( -C  +  /c2  +  4o2 ) ,  x  =  ±1 

(2.27)  XV  -  n1  -  p1  -  C  . 


In  the  following  chapters  we  will  carry  out  our  asymptotic  analysis  for  (2*19) 

-  (2.24)  for  X  ♦  0  which  corresponds  to  very  high  doping  (C  ♦  •).  This 

usually  is  the  case  for  modern  devices#  For  the  set  of  parameters  given  in 

—  17  -3 

Table  2  of  Chapter  1  and  a  maximal  doping  c  »  10  cm  X  would  take  the 

value  X  «  10  3.  We  now  use  the  fact  that  C(x)  is  an  odd  function  to 
sisplify  (2.19)  -  (2.24):  If  we  substitute 

(2.28)  t(x,t)  *  4(x,t)  -  j  (♦1  (t )  +  t^ft))  • 

Equations  (2.19)  -  (2.21)  remain  unchanged. 

The  boundary  conditions  (2.22)  become 

(2.29)  #(1,t)  -  -*(~1,t)  -  j  (^(t)  -  t_t(t))  =:  *t(t> 

and  the  initial  condition  for  in  (2.24)  has  to  be  replaced  by 

(2.30)  f(x,0)  -  ^(x)  ^(x)  -  j  (*.,(0)  +  ^(O))  . 


(2.19)  -  (2.23)  can  now  be  simplified  by  the  "Ansatz" 


m 


(2.31)  $(-x,t)  -  -♦(x#t),  n(-x,t)  =  p( x, t ) ,  Jn(-x,t)  *  Jp(x,t) 

For  the  rest  of  this  paper  we  assume  that  this  "Ansatz"  is  compatible  with  the 
initial  conditions  (2.24) •  Thus  we  make  the 

(2.32)  Assumption?  n1  and  p1  in  (2.24)  satisfy 


(2.33) 


♦x(x)  =  -♦I(-X)  +  *I(1)  +  «|»I(-1),  nX(x)  =  pX(-x) 


The  "Ansatz"  (2.31)  gives  the  conditions 


(2.34) 


♦(0,t)  =  0,  n(0,t)  *  p(0,t) ,  Jn(0,t)  =  Jp(0,t) 


Thus  we  obtain  the  simplified  problem  on  [0,1] 

<2.35)  -  P  -  C 

(2.36)  nt  -  *  R 

(2.37)  p.  =  -J  +  R 

”x 

(2.38)  Jn  -  nx  -  n*x 

(2.39)  Jp  -  -px  -  p*x 
x  e  [0,1] 

(2.40)  *(0,t>  =  o,  n(0,t)  *  p(0,t) ,  J  (0,c)  =  J  (0,t) 

n  p 

(2.41)  ♦dpt)  -  ^(t),  n(  1  ,t )  =  ~  (C(1)  +  /C(1)2  +  4a2)  , 

p(1,t)  -  j  (“CM)  +  /C(1)2  +  4a2) 

(2.42)  ♦(x,0)  *  ♦I(x)  ,  n(x,0 )  =  nX(x),  p(x,0)  =  pX(x) 

(where  for  simplicity  we  have  written  instead  of  i|>,  and  if1). 

Remark :  The  reason,  why  (2.35)  -  (2.42)  is  simpler  than  (2.19)  -  (2.24),  is 
that  the  solution  of  (2.19)  -  (2.24)  will  exhibit  an  internal  layer  at  x  =  0 
which  can  be  treated  as  a  boundary  layer  in  (2.35)  -  (2.42). 


v* *  ' 


-v  *  /.  V  w\ 


*  .  * 

*  O  -  * 


‘  /  A.  A 
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3.  Review  of  the  Singular  Perturbation  Analysis  for  the  Steady  State  Case 

We  now  review  the  singular  perturbation  analysis  carried  out  in  Markowich 
and  Ringhofer  (1982J  and  Markowich  et  al  [1983a]  for  the  steady  state  problem 

(3.1)  *2+xx  -  n  -  p  -  C 

(3.2)  Jn^  =»  R  ,  Jn  -  -  n+x 

(3.3)  JPx  =  "R  '  Jp  "  “Px  -  P*x 

(3.4)  R  *  (np  -  a2)(n+p  +  2a)”1 

(3.5)  *(0)  -  0  ,  *(1)  - 

(3.6)  n  =  j  (C  +  /c2  +  4a2),  p  -  j  (-C  +  /c2  +  4a2),  x  -  tl  . 

The  purpose  of  this  review  is  twofold:  On  one  hand  there  will  be  a  fair 
amounbt  of  analogy  between  the  asymptotic  behavior  of  (3e1)  -  (3.6)  and  the 
behavior  of  (2.35)  -  (2.42)  on  the  "slow"  time  scale.  On  the  other  hand  the 
initial  values  #  n 1  and  p1  in  (2,42)  will  usually  be  the  solutions  of 
the  steady  state  problem  (3.1)  -  (3.6)  and  it  is  thus  important  to  know  their 
structure.  Markowich  and  Ringhofer  [1982]  derived  an  asymptotic  expansion  for 
the  problem  (3.1)  -  (3.6)  in  powers  of  X.  This  expansion  is  of  the  form 

*  -  ^  i 

(3.7)  w(x,X)  ~  \  [wi(x)  +  w  (x/X)]X 

i*0 

(3.8)  w  *  (♦,n,p,Jn,Jp)T  . 

Thus  the  solution  of  (3.1)  *  (3.6)  can  be  represented  as  the  sum  of  a  smooth 

00 

part  (which  has  an  asymptotic  expansion  2.  w^(x)X  )  and  a  boundary  layer 

i=0 

term  at  x  =  0  which  varies  on  the  space  scale  x/X  and  has  an  asymptotic 

00 

*  i 

expansion  of  the  form  )  w  (x/X)X  .  The  solution  of  the  reduced  problem 

i=0 

Wg  satisfies  the  equations 


-10- 


0  “  "o  -  P0  '  c 

\  ■  *<VV  '  \  “  -*(V*o> 

X  X 

A#  <v  <V  y  A#  <w 

J  ”  n.  -  n.9.  ,  J  *  -p.  -p,*«  • 

n„  0  0  0  p_  0  *Tr0 

Ox  x  0  x  x 

(3.10)  -  (3.11)  is  a  system  of  2  second  order  ordinary  differential  equations 

(obtained  by  inserting  the  egressions  (3.11)  for  the  current  densities  JR 

and  Jp  into  the  continuity  equations  (3.10))  coupled  to  the  algebraic 

condition  (3.9).  Thus  we  can  impose  two  boundary  conditions  at  each  of  the 

boundaries  x  ■  0  and  x  *  1.  These  boundary  conditions  are 


(3.9) 

(3.10) 

(3.11) 


(3.12)  exp[-?Q(0)]n0(0)  -  exp[*0(0)]p0(0) 

(3.13)  J  (0)  -  J  (0) 

n0  P0 

(3.14)  ?0(1)  *  *1  '  Pi  "  J  (mC  *  +  *°2)  ‘ 

Note  that  the  boundary  conditions  (3.14)  together  with  (3e9)  imply  the  third 

1  2  2 

boundary  condition  a  m  j  (C  +  /c  +4o)  in  (2.41).  This  is  why  no  boundary 

A 

layer  occurs  at  x  *  1.  The  zeroth  order  layer  term  w^  satisfies  the 
relations 

(3.15)  nQ(0  -  n0(0){exp[f0(C)]  -  1},  P„(C)  “  p0(0){exp(-*0(c)l  -  1} 


Thus  for  ths  ssroth  order  approximation  of  w  we  have  a  boundary  layer  at 
x  ■  0  in  the  variables  ♦ ,  n  and  p.  Jn  and  Jp  do  not  exhibit  a  boundary 
layer  at  x  ■  0. 


.  /  V  V  V  V  '.*  V  *.* *.  *  ■ 

.  */  *.  *,•  *§.  •  •  *#* 


>totic  teglygif  for  the  Time  Dependent  Problem 


In  this  chapter  m  carry  out  our  asymptotic  analysis  for  the  solution  of 
(1.1)  -  (1.12)  as  X  tends  to  zero.  As  we  will  see  the  choice  of  initial 
conditions  crucially  influences  the  asymptotic  behavior  of  (1.1)  -  (1.12).  He 
assume  that  the  initial  values  f1?  n1  and  p1  have  an  asymptotic  expansion 
of  the  form 


(4.1) 


wZ(x,X)w  Y  CwJ(x)  +  w1^)]!1 
i-0  1  A 

w1  ”  (#1,nI#pI)T»  lwI(C)l  <  C^expt-CjC) 


This  choice  is  of  course  motivated  by  the  structure  of  the  solutions  of  the 
steady  state  problem  (1.17)  -  (1.20).  First  we  derive  our  asymptotic 
expansion  for  the  solution  of  (1.1)  -  (1.12)  which  varies  on  the  slow  time 
scale  t^  i.e. 


(4.2) 


w(x,t,X)  ~  Y  tw.(x,t)  +  v  (f'tHA1  . 
i-0 


In  order  for  this  expansion  to  be  valid,  it  will  be  necessary  to  impose 
certain  restrictions  on  the  initial  values  n1  and  p1.  These  conditions 

will  be  satisfied  if  »  n1  and  p1  are  the  solutions  of  the  steady  state 
problem  (1.17)  -  (1.20).  The  zeroth  order  term  wQ  of  the  outer  solution 
satisfies 


(4.3) 


0  *  n0  "  p0  "  C 


(4.4) 


(a)  nft  -  J  -  R(nn,pft),  (b)  prt  -  -J  -  R(nft,pn) 


°t  no 


°t  P0 


(4.5)  (a)  J  -  nQ  -  nQ*0  .  (b)  J  *  ~p0  -  pQ*Q  . 

XX  X  X 

Substracting  (4.4) (b)  from  (4.4) (a)  and  replacing  n  by  P0  +  C  gives 


(4.6) 


•  -  3L  .  3L  -  ar  ♦  J  -  C  -  <2pn  ♦  c)*r 


f  A  a  W  A  ^  V  ■  w  —  V 

0  0  n„  p„  x 

x  0  *0 


(4.7) 


-3n  -  R(p0  +  C,pQ),  Jp  -  -pfl  -  pft*„  . 


0  ‘'(To 

x  x 


(4.6)  is  now  a  second  order  elliptic  equation  for  ^  coupled  to  the 
parabolic  equation  (4.7)  for  pg«  Thus  we  can  impose  two  boundary  conditions 
at  each  boundary  x  ■  0  and  x  ■  1.  At  x  -  1  we  impose  the  conditions 
(4.8)  f0(1,t)  -  #t(t),  pQ(1,t)  -  j  (C  +  /C2  +  4a2)  . 

Mote  that  -  as  in  the  steady  state  case  -  (4.8)  together  with  (4.3)  implies 
the  third  boundary  condition  ng(1,t)  “  (C(1)  +  /C(1)^  ♦  ♦<»*)  in  (2.41). 
Thus  we  have  satisfied  all  three  boundary  conditions  at  x  *  1  with  Wq. 

This  is  the  reason  why  there  is  no  boundary  layer  at  x  *  1.  Inserting 
Wq(x)  +  wq(0  into  equations  (2.35)  -"(2.42)  gives 


(4.9) 


no  *  p0  'X 


(4.10) 


A[n0  +  R(n0  ♦  nQ.  p„  +  pQ)  -  R(JT0,P0)1 


(4.11) 


-  R(S0  +  »0.  pQ  +  p0)  +  R(^0.P0)1 


(4.12) 


(4.13) 


-  ”ot  -  <».  * 

«  -  -5,  -  «. ♦  ■  **  • 


X  ♦  0  gives 


(4.14)  «  n0  -  p0 


(4.15)  J  “  J 

n0  P0_ 


(4.16)  nQ  -  [nQ  +  nfl(0,t)U>0 


P0  *  <P0(0't)  +  "  0 


14- 


Since  «v0(C,t)l  <  C1(t)expl-C2(t)Cl  has  to  hold,  w(«*,t)  -  0  holds.  Thus 
we  can  integrate  (4.15)  and  (4.16)  and  obtain 


(4.17)  J  =  J  SO 

n„  p- 
0  ^0 

(4»18>  nQ(C,t)  -  nQ(0,t){exp(t0(C«t)]  -  1> 

(4.19)  pQ(C,t)  -  P0(0,t){exp[-i0(c,t)]  -  1} 

(4.20)  ♦„  -  n0  -  pQ  . 

Inserting  Wq  +  w^  in  the  boundary  conditions  gives 

(4.21)  V°»t)  “  -?0<0,t) 

(4.22)  n0(0,t)exp[f0(0,t)l  -  p0(0,t)exp[-J0(0,t) ] 

(4.23)  3f  (0,t)  -  3F  (0,t)  . 

n0  P0 

Coabining  (4.21)  and  (4.22)  gives 

(4.24)  n(0,t)exp[-f(0,t)]  ■  p(0,t)exptf(0,t)] 

which  together  with  (4.23)  gives  the  boundary  conditions  for  the  reduced 
problem  at  x  •  0.  We  now  investigate  the  initial  conditions  for  the  reduced 
equations  (4.3)  -  (4.5).  Since  equation  (4.6)  is  elliptic  it  has  to  hold  also 
at  t  ■  0.  This  gives  the  condition 

(4.25)  [Cx  -  (2?J  ♦  C)?Jlx  -  0 

(4.26)  nj(0)expt-^j(0))  -  pj(0)exp[?j(0)] 


(4.27)  ?J(1)  -  ^,(0) 

on  the  seroth  order  term  of  w  (x,X).  Also  the  layer  equations  (4.14)  - 
(4.16)  have  to  hold  for  t  *  0  since  they  are  time-independent.  This  gives 
the  condition 
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*S 

(4.26) 

aT 

nj(t>  - 

•» 

A  — 

(4.27) 

1 

u 

-1  o 
Cm 

M 

r, 

il 


>u 


n?(C)  -  nJ(0){exp[+*E(C)J  -  1} 


pJ(C)  -  p*(0){exp[-*J(C)]  -1> 


*1 

on  Wq.  Equation  (4.14)  has  to  be  satisfied  because  of  the  compatibility 


condition  (2.28).  Using  the  theory  outlined  in  the  previous  chapter  we  see 
that  conditions  (4.24)  -  (4.27)  are  satisfied  if  w1  is  the  solution  of  the 
steady-state  problem  (3.1)  -  (3.6).  Thus,  if  conditions  (4.24)  -  (4.27)  are 


~I  A1 

satisfied  we  have  determined  the  zeroth  order  term  wQ  +  wQ  of  the  expansion 


(4.2).  We  summarize  this  in  the 
(4.28)  Theorem : 


If  wJ(x,X)  has  an  asymptotic  expansion  w1  —  ^(w*(x)  +  w*(^-)  )X*» 

i»0  1  1  A 

■  (#I,nI,pI)T  and  wj  and  w*  satisfies  the  conditions  (4.24)  -  (4.27) 


(as  it  is  the  case  if  w1  is  the  solution  of  the  steady  state  problem  (4.1) 

»  M  A  «  1  _ .  •  _  _ .  •  « . _  .A  #  »  \  .  I  •  ■  .  •  .  . 


(4.6))  then  the  zeroth  order  term  w^(x,t)  satisfies  the  equations 

(4.29)  °  -  ?0  ,  ?0  -  3  +  3  «  cx  -  (2p0  +  C)^0 

x  n0  p0  x 

(4.30)  5.  -  -3  *  k<5.  ♦  c,p„>, 


;o  -  -Jp„  *  wfo  ♦  c'?o>'  \ 
t  o 

X 


<V  a#  y 

■Pq  ~  p<ro 


together  with  the  boundary  conditions 

(4.31)  [p(0,t)  +  C(0)]exp[-5»Q(0,t)]  -  pQ(0, t)exp ($0( 0,t) ] 


(4.32) 


JQ(0,t)  *  2Jp  (0,t) 


(4.33)  ^f0d,t)  -  ^(t),  pQd,t)  -  ~  (-C(I)  +  /C(1)2  +  ♦a2) 


and  the  initial  conditions 
(4.34) 


P0(x,°)  »  pQ(x) 


The  zeroth  order  layer  term  w^(C,t)  satisfies  the  relations 
(4.35)  nQ(C,t)  -  nQ(0,t) [expC4Q(C,t)]  -  1] 


uU 


(4.36) 


(4.37) 

(4.38) 


pQ(C#t)  -  p0(0,t) [exp[-f0(C,t)]  -  1] 

A  A 

J  -  J  -  0 
n.  p. 

0  0 

♦0  ■  n0  ”  P0#  V0'0  ■  -^0(o,t),  t.(»,t)  »  0  . 

w  o  u 


So  we  have  determined  the  zeroth  order  approximation  to  the  solution  of  (2.35) 

-  (2.42)  in  the  case  that  the  initial  functions  ♦*#  n1  and  p1  are 

solutions  of  the  steady  state  problem.  This  expansion  varies  in  time  on  the 

timescale  t.  Thus  the  reaction  time  of  the  p  -  n  junction  (the  time  it 

takes  the  system  to  reach  another  steady  state  if  we  vary  the  bias)  will  be  of 

this  order  of  magnitude.  Bssential  for  the  derivation  of  this  asymptotic 

solution  was  the  fact  that  (4.24)  -  (4.27)  was  satisfied  for  the  initial 

values.  Although  this  will  be  the  case  for  all  practical  purposes  (since  we 

always  start  with  a  steady  state  solution)  it  is  of  some  interest  to  see  what 

happens  if  conditions  (4.24)  -  (4.27)  are  violated.  This  is  the  case  if  we 

solve  (2.35)  -  (2.42)  numerically  and  our  numerical  solution  is  "polluted"  by 

a  discretization  error.  For  this  purpose  it  is  convenient  to  rewrite  our 

problem  (2.35)  -  (2.42) s  Me  elliminate  n  from  (2.36)  -  (2.38)  by 

2 

substituting  n  *  X  f  +  p  +  C  in  (2.36)  -  (2.38)  and  subtracting  (2.37) 
from  (2.36).  This  gives 


(4.39)  X2^  *  »x.  J  -  Jn  +  *  C,  ♦  X2^ 


(4.40)  pfc  -  -J  ~  R,  J  m  “Px  “  p#x 


We  now  supplement  our  expansion  (4.2  )  by  a  term 
the  fast  time  scale  T  -  tX"*,  Y  >  0  (where  y 


(2p  +  C  +  X2t  )♦ 

r  TXX  X 


w(t,x,X)  which  varies  on 
still  remains  to  be 


determined ) 


(4.41) 

(4.42) 


“  ~  -  i 

w(x,t,X)  ~  \  [w^X/t)  +  +  W^X/TJJX 


i-0 


lwi(x/T)l  <  C1 (x)exp(-C2(x)T)/  i  -  1/2/ 


5.  Flnlf  Difference  Methods 

In  this  chapter  we  develop  a  finite  difference  scheme  for  the  solution  of 
problem  (1.1)  -  (1.12).  Severe  stability  problems  arise  because  of  the 
presence  of  the  fast  time  scale  (t/X  )  in  (4.47).  To  illustrate  this  we 
first  consider  a  "naive"  discretisation  of  (2.35)  -  (2.42):  First  we  only 
dlscretise  (2.35)  -  (2.42)  in  time  and  leave  the  space  variable  x 
continuous.  Thus  we  define  a  sequence 

(5.1)  {tfc,  0  -  t0  <....<  tfc  <....  },  :«  tk+1  -  tx,  k  -  0,1,... 

and  denote  the  function  v(x, )  with 

(5.2)  vk(x)  ■  v(t.  ,x),  v  -  (#,n,p,J  ,J  )T,  k  -  0,1,...  . 

jc  n  p 

As  a  first  approach  one  would  perhaps  try  to  solve  equations  (2.35)  -  (2.42) 
by  a  Crank -Nicholson  scheme  and  solve  Poisson's  equation  (2.35)  implicitly  on 
the  new  time  level.  This  would  give 


(5.3) 

x  V*'  -  1 

Txx 

„kt1  -  pk*'  -  c 

(5.4) 

k+1  k 
n  -n 

Atk 

5  (C  *  <  > 

2  “x  “x 

(5.5) 

k+1  k 

P  "P  - 

Atk 

-  i  (Jkk’  +  3*  ) 

2  V  D  D  * 

*x 

(5.6) 

J4  -  n4  - 
n  x 

n4** 

"  ♦« 

i  ■  k,  k  +  1 

(5.7) 

£  £ 

J*  *  -p* 
p  *x 

.  JL.% 

-  P  *x  • 

As  we  saw  in  chapter  4  the  part  of  the  potential  ♦  varying  on  the  fast  time 
2 

scale  (t/X  )  satisfies  the  equation 

■  Icx  •  <"4*>*x1x  • 

k+1  v 

Time  differencing  equation  (5.3)  and  inserting  the  expressions  n  -  n  , 
pk+1  -  pl*  froB  (5.4)  -  (5.7)  yields 
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i 


J 


£ 


* 

.1 


*5 


•  k+1  _  .k 

(5.8*2  -£L_ - —  »  \  [2C  -  (nk+1  +  pk+1)*k+1  -  (nk  +  pk)*k]  +  0(X2) 

2  X  X  XX 


Thus  away  from  the  boundary  layer  at  x  =  0,  <|>xx  will  behave  like  the 


|» 

solution  y  of 


(5.9) 


%2  yk+1  -  yk  1  ,  k  k  .  k+1  k+1,  ^k 

x  h^1— !<•>  **  »  >*'•’*>  0 


it. 


which  oscillates  wildly  unless  ~  <  const  holds.  To  derive  a  scheme  which 

X2 


X  :*  {x.  :x.  ■  ih,  i  *  0(1)N}#  h  *  ^  . 

11  n 


does  not  need  this  severe  restriction  on  the  time  steps  we  rewrite  (2.35)  - 
(2.42):  He  proceed  as  follows:  He  first  discretize  (2.35)  -  (2.42)  in  the 
space  variable  x.  Thus  we  define  a  mesh  X  and  a  meshsize  h  by 
(5.10) 

(For  the  sake  of  siag>licity  we  assume  a  uniform  grid.  The  generalization  of 
our  results  to  nonuniform  grids  is  straightforward. )  Realistically  we  have  to 
assume  that  h  »  X  holds. 

T  t  i  T 

He  approximate  (4»n,p)  at  xA  by  the  gridfunction  • 

i  »  0(1 )N>  He  approximate  the  currents  (Jn*J_)  at  x  i  !“  (xj_  +  xi+i^2 


i+ 


by  {(J  ,  J  ),  i  ■  0(1)N-1}.  For  a  gridfunction  {z.  *  i  *0(1 )N)  we 

n  ^  d  i 


1 


i*  2  a 


define  the  discrete  operators  D+#  D.  and  M  by 
(5.11)  m  *1+1  "  *i #  **  "  *i—i 


1 


(Mz)A  -  ?  (sA  +  zi+1)  . 


(5.12) 

He  discretize  (1.1)*  (1.2)  and  (1.3)  by 


(5.13) 


2.-2 


X  h 


<D+D>>i  "  "j.  -  Pi  “  C(xA) 


(5.14) 


dt  ni  ”  -Jn 


-  J 


n 


)A 
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(5.15) 


ft  “i  '  -lJp  .  -  Jp  )^>  ' 

1+  5  1‘  2 

For  simplicity  we  set  the  recombination  rate  R  «  0.  For  the  discretization 
of  the  current  relations  (2.38),  (2.39)  we  use  an  exponentially  fitted 
scheme.  This  scheme  was  first  proposed  by  Scharfetter  and  Gummel  [1969]  and 
analyzed  in  a  singular  perturbation  framework  by  Markowich  et  al  [1983b]: 


(5.16) 

(5.17) 


J 


0(z)  : 


'  [0(D+*i)(D+ni)  -  (Mni)(D+<|>1)]/hi 
•  I-CtD^  )(D+p1  )  -  (Mpi)(D+t|>i)] 

^  coth  (z)  . 


The  boundary  conditions  transform  to 


(5.18)  *0(t>  -  0,  nptt)  -  P0(t),  Jn  (t)  -  Jp  (t) 

2  2 

(5.19)  tN(t)  -  ^(t),  p^(t)  -  ^  (-C(I)  +  /c(1)2  +  4a2),  i^(t)  =  pN(t)  +  C(1) 
and  as  initial  conditions  we  impose 

(5.20)  ^(0)  -  ni(0)  *  nj,  Pi(0)  -  pj,  i  =  1(1  )N~1  . 

Where  for  the  sake  of  compatibility  we  assume 

(5.21)  X2h"2D+D_*J  *  ni  “  pJ  “  C(x  ),  i  *  2(1)N-1  . 

(5.13)  -  (5,20)  is  an  initial  value  problem  for  a  system  of  2(N-1)  ordinary 
differential  equations  coupled  to  the  N  -  1  algebraic  equations  (5.13).  We 
now  eliminate  the  algebraic  equations  (5.13)  by  differentiating  (5.13)  with 
respect  to  t  and  inserting  for  nt  and  pt  from  (5.14)  and  (5.15).  We 
obtain 
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(5.22) 


S 

IS 


X“h  *D 


'+D-  It  '  (J..  i  *  J.  ,V> 


i+  -  i-  - 
2  2 


(5.23) 


J.+  1  "  Jn  ,  +  Jp  “  [0(Mi,(D+ni  -  D+Pi)  - 

1  i+  2  i+2 


(Mn^  +  Mp^tD^H/h 


We  now  eliminate  n^  from  (5.23)  by  using  (5.13)  and  the  boundary  conditions 


(5.18),  (5.19) 


(5.24) 


n^  ■  pA  +  v^Cix^)  +  w^X^h”^D+D_^ 


(5.25) 


r  1  i  f  0,N  f  1  i  +  0 

*i  *■  1  vi  ■  1 

^  0  i  -  0,N'  ^  0  i  -  0  . 


Thus  we  obtain  the  system 


(5.26)  XV(D+D_  fj.  ♦1)i  -  tj  ,  -  J  ,) 


i+  -  i-  - 
2  2 


(5.27)  hJ  t  -  {o(D+*i)D+(C(xi))  -  M(2p+C)1 (D+*A )±  + 


i+  2 


+  X2h"2[0(D+*i)D+(wD+D_t(>)1  -  (M(wD+0_*))i(D+*)1]} 


which  approximates  the  fourth  order  differential  equation 


(5.28) 


X  *xxt  -  Jx 


(5.29) 


J  -  C  -  (2p+c)*  +  X2*  -  X2*  * 

X  X  XXX  XX  X 


The  problem  consisting  of  (5.15),  (5.17),  (5.26),  (5.27)  together  with  the 


initial  conditions  (5.20)  for  pA  and  is  equivalent  to  the  original 

problem  (5.13)  -  (5.17)  if  we  define  n^  by  (5.24). 


Inserting  (5.24)  into  (5.23)  we  obtain 


(5.30)  X2  (D+D_  f)  -  D_(o(D+*)(D+C)  -  M(  2p+C)  (D+*) ] ±  +  O(^) 

h 


(5.31) 


p)i  »  h  D_[0<D+*HD+p>  +  (Mp)(D+*)]i 


To  avoid  the  stability  problems  outlined  at  the  beginning  of  this  chapter 


(5.30)  has  to  be  discretized  by  backward  differencing.  (5.31)  can  be 


discretized  by  the  trapezoidal  rule  (or  the  Crank-Nicholson  scheme).  The 


terms  abreviated  with  0(A  h  )  in  (5.30)  can  be  taken  at  the  previous  time 


step  since  we  assumed  Ah  <<  1.  Thus  we  approximate  n^,  p^,  J  ^ , 


*  ^  ^  ^  ^  a  k  k  k  k  k 

J  at  t  *  tv  by  ♦ . ,  n  .  p  #  J  J 

P1  iiii+iP1 

i+2  2  i+I 


and  obtain 


(5.32) 


X2(h,At)_1(D  D  (*k+1  -  *k))  -  Jk+1  -  Jk+1 

1+  2  A-  2 


(5.33) 


,A.I-1/  k+1  k*  1  /  Tk+1  _k+1  ,  _k  ,  _k  \ 

(4t)  (Pi  -  Pt)  -  j  (Jp  -Jp  *  Jp  ♦  Jp  )/h 

I  i_  2  i+  5  l-  2 


V  k+1  V+1  k+1 

hJK  1  -  <J(D+**  )(D+C)i  -  M(2p*  +  C)i(D+*K  ,)i  + 

i+  2 


(5.34) 


+  X2h"2  to(D  *k)DA(w(D  D  *k)).  -  M(w(D  D  i|»k) )  .  (D.ipk)  . 

+  i  +  +  —  i  +  —  i  +  x 


(5.35) 


Jk  -  -0(D+tJ)(D+pk)i  -  (Mpk)i(D+<|.k)i 

i+2 


(5.36) 


.k+1  .  .k+1  ,  tJ_  . 

♦o  "  °'  *N  *Vtk+1) 


(5.37)  2Jk+1-  Jk+1  pjj+1  -  Pb(tk+1)  =  \  t-C(1)  +  /C(1)2  +  4a2]  . 


This  corresponds  to  solving  the  differential  equations 


(5.38)  7-  (♦k+1  -  )  -  [C  -  ( 2pk+ 1  +  C)*k+1  +  X2||»k  -  X2*k  i|>k] 

At  XX  XX  X  X  XXX  XX  X  X 


'  V  *  ■ 


(5.39) 


IA«.l“1fJC+1  _  *  1  f„k+1  +  r,k+1*k+1i  _  1  Jk 

(At)  (p  -  P  >  *  2  (px  +  P  *x  }x  2  Jp 


at  each  time  step,  which  represent  a  coupled  system  of  two  nonlinear  equations 
for  f  and  p.  (5.33),  (5.35)  (or  (5.37))  represent  a  second  order  time 
discretization  for  (5.15),  (5.17)  whereas  (5.32),  (5.34)  only  is  a  first  order 
time  discretization  for  (5.22),  (5.23). 

2 

The  time  derivative  in  (5.22)  is  multiplied  by  X  .  Thus  the  method  is 
actually  convergent  of  order 

0[(At)2  +  A2At] 


where  At  denotes  the  maximal  time  step.  Since  for  practical  purposes 
At  >>  X  will  hold  we  can  say  that  the  method  is  actually  of  second  order  in 
time.  A  detailed  analysis  of  the  stability  and  convergence  properties  will  be 


the  topic  of  a  subsequent  paper 


6.  Numerical  results 


In  this  chapter  we  present  numerical  results  for  two  test  problems  which 
were  solved  by  the  technique  outlined  in  the  previous  chapter.  The 
confutation  was  performed  on  a  VAX-780  computer  in  double  precision  (which 
provides  14  significant  digits).  The  physical  parameters  were  taken  from 

Table  2  in  Chapter  1.  We  took  a  constant  doping  profile  on  the  interval 

—  17  —3 

[0,£].  The  maximal  doping  (C  in  Chapter  2)  was  taken  to  be  10  cm  .  The 

characteristic  length  21  of  the  device  was  taken  to  be  5y.  This 

—3 

corresponds  to  X  *  10  •  As  initial  values  the  solution  of  the  equilibrium 

problem  (that  is  the  steady  state  problem  with  zero  bias  applied)  have  been 
taken  on  the  interval  [0,1].  21  equally  spaced  gridpoints  were  taken  in 

the  x  direction  (h  *  0.05).  As  time  steps  we  used  At^  =  At  =  0.1, 
k  *  0,1,2,...  •  (No  stability  problems  were  encountered  when  we  varied  At.) 
Two  cases  have  been  simulated:  First  we  switched  from  the  equilibrium  to  a 
reverse  bias  of  0.25V.  In  the  other  case  we  switched  from  0V*  (equilibrium) 
to  a  voltage  of  0.25V  forward  bias.  This  corresponds  to  a  variation  of  the 
boundary  condition  ♦^(t)  in  (2.41)  from  ^  *  16.12  (~  ff.V)  to  ^  -  21.12 
(~  0.25V  reverse  bias)  in  figures  1,  2  and  3  an  to  a  switch  from  tp^  =  16.12 
(~  0.V)  to  *■  11,12  (~  0.25V  forward  bias)  in  figures  4,  5  and  6. 

We  plotted  iD  and  the  hole  density  p  on  the  interval  [0,1].  The 
electron  density  n  would  be  given  by 

n(x, t)  -  p(x,t)  +  C(x)  +  X21>xx(x,t) 

Thus  away  from  x  *  0  it  differs  visually  from  p  only  by  the  (constant) 

function  C(x).  To  obtain  a  picture  of  the  full  solution  on  [-1,1]  ip,  n 

and  p  would  have  to  be  continued  according  to  <p(-x)  *  -ip(x),  n(-x)  ■  p(x). 

In  figures  3  and  6  we  plotted  the  value  of  the  total  current 

J(t)  -I1,  J  (x,t)  +  J  (x,t)dx 
'-in  p 

-25- 


H 

j 

.>4 


as  a  function  of  time. 

In  detail  tha  figures  1  -  6  shows 
Figures  1  -  3s  Switch  from  ?.V  to  0.25V  forward  bias 

Figure  Is  f(x,t),  Figure  2s  log1Q (p(x,t) )  , 

Figure  3s  log^Q(J(t)) 

Figures  4  -  6s  Switch  from  ff.V  to  9.25V  forward  bias 

Figure  4s  #(x,t),  Figure  5s  log1Q (p(x,t) ,  Figure  6s  log1Q(J(t)) 
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